Commensurate lock-in and incommensurate supersolid phases of hardcore bosons on 

anisotropic triangular lattices 
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We investigate the interplay between commensurate lock-in and incommensurate supersolid phases 
of the hardcore bosons at half-filling with anisotropic nearest-neighbor hopping and repulsive inter- 
actions on triangular lattice. We use numerical quantum and variational Monte Carlo as well as 
analytical Schwinger boson mean-field analysis to establish the ground states and phase diagram. 
It is shown that, for finite size systems, there exist a series of jumps between different supersolid 
phases as the anisotropy parameter is changed. The density ordering wavevectors are locked to 
commensurate values and jump between adjacent supersolids. In the thermodynamic limit, how- 
ever, the magnitude of these jumps vanishes leading to a continuous set of novel incommensurate 
supersoild phases. 
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Introduction. — Supersolid phases with coexisting long- 
range diagonal and off-diagonal orders have long been 
recognized as interesting conceivable ground states of su- 
perfluid systems in the presence of mobile vacancies 
Recent experiments have also reported possible evidence 
for such a state in 4 He [2]. Another, relatively new, 
route to supersolid phases has recently been explored in 
the hardcore boson models on frustrated two-dimensional 
(2D) lattices [3|, In such lattices, supersolid phases 
may arise due to intricate competition between kinetic 
and interaction energies, and as a result present an ex- 
cellent playground for discovery of possible novel univer- 
sality classes of quantum phase transitions [B|, 0] . Further 
interest in these systems stems from potential realization 
of these models in ultracold atomic systems on optical 
lattices Q- 

The diagonal (density) orders of the supersolid phases 
discovered so far have been restricted to commensurate 
orders. For example, a commensurate supersolid was dis- 
covered for the hardcore bosons on the isotropic trian- 
gular lattice Q with the density ordering wave vector 
Qo = (47r/3, 0) along with non-vanishing superfluid or- 
der parameter A qualitative understanding of such 
a supersolid state can be obtained by considering the 
deviation from a commensurate boson filling fraction at 
which the ground state is a perfect Mott crystal with a 
long-range diagonal order. The additional particles or 
holes resulting from such a deviation condense to pro- 
duce a superfluid while retaining the backbone of the 
existing Mott solid. This leads to the coexistence of off- 
diagonal (superfluid) and diagonal (Mott) orders. One 
may expect that incommensurate version of the super- 
solid phases may arise when there are more than one 



competing interactions or length scales in the system. If 
such a phase exists in lattice models, this may be a much 
closer analog of the supersolid phases originally proposed 
for the continuum. 

In this letter, we investigate possible presence of in- 
commensurate supersolid phases of the hardcore bosons 
at half-filling with anisotropic nearest-neighbor hopping 
and repulsive interactions on triangular lattices. We es- 
tablish the ground states and phase diagram using several 
different and complementary methods, namely numerical 
quantum and variational Monte Carlo (QMC) techniques 
as well as analytical Schwingcr-Boson mean-field theory. 
It is found that, for finite size systems, there exist a series 
of jumps between different supersolid phases as a func- 
tion of the anisotropy parameter. The (density) ordering 
wavevectors of these phases are pinned to commensurate 
values and jump upon entering a nearby supersolid phase. 
The ordering wavevectors assume every single commen- 
surate values for an interval of the anisotropy parameters. 
In the thermodynamic limit, however, we find that these 
ordering wavevectors become a continuous function of 
the anisotropy parameter, leading to a continuous set of 
incommensurate supersolids. We emphasize that such a 
continuous set of supersolid orders represent exciting dis- 
covery of novel quantum structures that have not been 
seen in previous studies of lattice boson systems. 

Lattice Model. — We begin with the following hardcore 
boson model on an anisotropic triangular lattice. 



H h = 



+ H - c -) + v a n i n i 



(i) 



where bi denotes the boson annihilation operator at site 
i and (ij) runs over the nearest-neighbor sites. The 
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FIG. 2: The equal time structure factor S(Q(n)) at the or- 
dering wavevector Q(n) as a function of n for different system 
sizes and temperatures. All other parameters are the same as 
in Fig. [T] Judging from the size dependence of the data, the 
structure factor is clearly finite in the thermodynamic limit. 



FIG. 1: The superfluid density along a? as a function of the 
anisotropy parameter n for different system sizes and temper- 
atures with Vi/ti = 10. The inset shows a similar plot for the 
superfluid density along aj.. Lines are guides to the eye. The 
number of steps for L — 12, 24, 48, 96 are 3, 5, 9, 17, respec- 
tively. The number of jumps between adjacent steps is hence 
2, 4, 8, 16, respectively. 



hopping tij and repulsive interaction Vij are given by 
hj — h(Vij =Vi) and ty = t 2 (V~ij = V 2 ) for the nearest- 
neighbor sites along the diagonal and horizontal bonds of 
a triangular lattice (if it is viewed as a square lattice with 
one additional diagonal bond per each plaquette). Here 
we shall fix the anisotropy parameter r\ = tijt\ = V2/V1. 
For 77 = 1, the model reduces to the well-known isotropic- 
triangular-lattice model. Such a hardcore boson model 
is also equivalent to an anisotropic spin-1/2 XXZ model 
via the well-known Holstein-Primakoff mapping 0] . 

In the isotropic case, the classical limit of this model 
(ii,2 = 0) has an extensive ground state degeneracy and 
power-law density-density correlations (or S z -Sz corre- 
lator in the XXZ model) at zero temperature [8]. The 
ground state degeneracy at the isotropic point is com- 
pletely lifted for rj < 1 and the system orders at Qi = 
(jr, 7r). On the contrary, the degeneracy is only partially 
lifted for rj > 1. Here each diagonal chain is ordered an- 
tiferromagnetically at Q2 = (tt, 0) but the chains can be 
shifted with respect to each other giving rise to 2 L ground 
states, where L is the linear system size. We expect the 
first (second) type of ordering in the quantum model for 
77 <C (3>)1 and call these phases solid I(II) for future 
reference. For the quantum model, when t% t 2 is turned 
on, it is well-known that the system exhibits a supersolid 
phase at 77 = 1 and large enough values of V\ jt\ H ■ The 
key point which we want to address in this paper is the 
fate of the supersolid phase when 77 7^ 1. 

Quantum Monte Carlo. — To address this problem, we 
perform QMC simulations using a plaquette generaliza- 



tion [§] of the Stochastic Series expansion (SSE) algo- 
rithm [1 01 ] . where the elementary lattice unit is a tri- 
angle; this results in improved efficiency for large val- 
ues of V\/t\. We measure the superfluid density along 
the diagonal (a2) and horizontal (ai) lattice directions 
by measuring the corresponding winding numbers W£. 



W£ /fiti, where is the inverse tem- 



Psi(2) 

perature. We also measure the equal time density- 
density correlator S(q)/N — (p qT p qT ), where p qT — 
{1/N) J2i Pit exp(iq • r^) and p iT is the boson density at 
site i and imaginary time r. 

We begin with the case 77 > 1. A plot of the super- 
fluid density p s i along the a 2 lattice direction is shown 
in Fig. [T] as function of 77 for different system sizes and 
temperatures. Notice that it exhibits a staircase struc- 
ture with the number of steps proportional to the system 
size. Also, when the system size is doubled, the number 
of jumps between adjacent steps doubles and the gap be- 
tween them decreases by half. The corresponding plot for 
p s i is shown in the inset of Fig. [T] As can be clearly seen 
from both plots, the superfluid density undergoes sev- 
eral discontinuous jumps before reaching zero at 77 ~ 1.8. 
Also, as shown in Fig.O the density-density correlator (in 
the thermodynamic limit) in each segment of 77 is finite 
at some ordering wavevector Q{rf) in the corresponding 
parameter range. Further, the ordering wavevector Q(tj) 
is a constant along the "plateau" and changes discontin- 
uously upon entering the next phase. Thus, for these 
finite size systems, the "plateaus" for 1 < 77 < 1.8 corre- 
spond to distinct supersolid phases with sharp transitions 
between them. 

Next we address the anisotropy dependence of the or- 
dering wavevectors. The evolution of Q = (k x ,k y ) as 
a function of 77 is shown in Fig. [3J Here we have cho- 
sen (k x ,k y /y3) = (7iibi + 7i2b2)/27rL, where bi = 
27T(l,l/ v / 3),b 2 = 2tt(0,2/V3) are the reciprocal lattice 
vectors and ni,n2 are integers. Comparing Fig. [3] and [1] 
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FIG. 3: The ordering wavevectors Q(r]) = (k x 
tion of n for different system sizes and temperatures. All the 
parameters are the same as in Fig.[T] Note that for odd k x L, 
k v cannot be zero and hence is off the axis. The data for 
L = 96 for r\ > 1.6 are not fully equilibrated, leading to the 
noisy behavior. 
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FIG. 4: Ground state phase diagram obtained by VMC as 
a function of n and Vi/ti for L — 24. The supersolid phase 
exists for Vi > 10ii and between 0.8 < r\ < 1.5. The inset 
shows the superfluid order parameter and the density-density 
correlators S(n,n) and S(2tt/3, 2n/3) as a function of 77 for 
Vi/ti — 9. The transition from the Mott solid I to supersolid 
phase occurs at n ~ 0.8. 



wc find that the system locks at rational wavevectors that 
are commensurate with the lattice in the "plateau" re- 
gions. The ordering wavevector goes from Q = (47r/3,0) 
at the isotropic point (77 = 1) to Qi = (tt, 0) at the tran- 
sition point to the solid II phase (77 ~ 1.8), picking all 
possible commensurate values in between. The discontin- 
uous jumps of Q(f?) between these commensurate values 
decrease in magnitude with increasing system size. The 
nature of the phase diagram for 77 < 1, with 0.84 < 77 < 1, 
turns out to be qualitatively similar. We again find a 
series of supersolid phase with the ordering wavevector 
pinned to commensurate values along the "plateaus" , be- 
fore the system reaches the Mott phase solid I at 77 = 0.84. 
The only difference comes from the fact that for 77 < 1, 
k x L is always even, so that k y remains pinned to zero 
throughout the phase diagram. 

Variational Monte Carlo. — We now supplement the 
QMC results with the VMC studies of the global phase 
diagram. Following Ref. we use a variational wave 
function: 

I*) = e -^*.i Vi >* nin *\$ ), (2) 
where |$ ) = (&L0H ) isthe 

non-interacting superfluid 
wave function and N is the total number of bosons. The 
components of the Jastrow potential, u^j = v(\Ri — Rj\) 
are independently optimized to take into account the cor- 
relations between particles at different sites. The varia- 
tional ground state energy decreases with larger number 
of Vij. For the present study, we have incorporated 15 
Vij parameters to make our result qualitatively and semi- 



quantitatively consistent with the QMC results. Stan- 
dard Metropolis algorithm is used to calculate the varia- 
tional ener gy. The method of statistical reconfiguration 
by Sorella [12j is employed to obtain the optimized pa- 
rameters. With the optimized wave function, we compute 
the superfluid density p(k) = eiiRi ~ Rj) ' k ( b t b j)i the 
density-density correlator 5(q), and obtain the ground 
state phase diagram from these quantities [l3l]. 

The VMC results are summarized in Fig. 2] where the 
ground state phase diagram is shown as a function of 77 
and Vi/fi for L — 24. Notice that a supersolid phase 
exists in the range 0.8 < 77 < 1.5 for Vi/ti ~ 8, which 
is qualitatively consistent with the QMC results. We 
have found that the upper limit of the phase boundary 
r] u ~ 1.5 depends on the system size, and progresses to- 
wards larger values with increasing L. The lower phase 
boundary 77; ~ 0.8 is virtually independent of the sys- 
tem size. For 77 < 77/ (> 77^), the system enters the Mott 
phase I(II) provided V\/ti > 8. For weaker interactions, 
the superfluid phase prevails for all 77. The inset of Fig. 
[4] shows the transition from the Mott I to the supersolid 
phase at 77 ~ 0.8 for L — 12 and Vi/ti — 9. Note that the 
density-density correlator S(2n/3, 2n/3) and the super- 
fluid order parameter p(k = 0) rise sharply while S(ir, tt) 
drops to zero around this point, signifying a transition 
from the Mott I to the supersolid phase. 

Schwinger Boson analysis. — Next, to obtain an ana- 
lytical understanding of the 77 dependence of the order- 
ing wavevector, and to determine the fate of the order- 
ing wavevector plateaus in the thermodynamic limit, we 
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FIG. 5: (Color online) Plot of q (ordering wavevector is given 
by Q — (q, q)) as a function of n for several finite system sizes 
(colored symbols) and in the thermodynamic limit (dashed 
line) as obtained from Schwinger boson mean-field analysis. 



carry out a mean-field Schwinger boson analysis of Eq. 
[TJ To this end, we start from the XXZ spin model de- 
scription of Eq. Q] [3] and rewrite these spins in terms 
Schwinger bosons: — a|cj and Sf = (ajai — c\ci)/2. 



Following standard procedure [lj] , we decouple the resul- 
tant Hamiltonian using the mean-field order parameters 
A a = {a\a i+a + c\c i+a ) and B a = (aiC i+a - Cia i+a ), 
where a — x,y for the horizontal/ vertical and a — z 
for the diagonal bonds emanating from site i, to get the 
Schwinger boson mean-field free energy 
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Here A k = A x cos(k x ) + A y cos(k y ) + A z cos(k x + k y ), 
B k = B x sin(fc x ) + BySm(k y ) + B z sin(k x + k y ), uj k = 
^/|A — Ak| 2 — B'£ is the spinon dispersion, and the pa- 
rameter A is used to enforce the constraint of 2S — 
a\oi + c\ci at the mean- field level. 

We obtain the values of mean-field variables at ground 
state by solving the saddle-point equations qX* = 
d Qg F = gf^ = 0. The minima of the spinon disper- 
sion u)(k) at k m i n — ±(q/2, q/2) gives the spin order of 
the XXZ model with ordering wave-vector Q = (q, q) 
as shown in Fig. 5, both at finite sizes and in the 
thermodynamic limit. In the thermodynamic limit, for 
small ti/Vi ~ 0.04, we find the supersolid phase for 
0.8 < 7] < 1.3 with continually varying ordering wave- 
vector Q(r/) = (q, q) as shown in Fig. [5j The correspond- 
ing spinon dispersion is gapless around both k — (which 
is a signature of superfluidity) and k = ±(q/2,q/2) 
(which signifies the long-ranged solid order). For finite- 
size systems, we find that the spinon dispersion at k = 
±(q/2, q/2) acquires a gap which decreases with increas- 
ing system size and vanishes in the thermodynamic limit. 



Such a gap of the spinon dispersion leads to the staircase 
behavior of the ordering wave-vector as shown in Fig. 
[5] Our analytical Schwinger boson results are in quali- 
tative agreement with both QMC and VMC results for 
finite system sizes and we therefore expect it to predict 
the correct behavior of the ordering wave-vector in the 
thermodynamic limit. 

To conclude, we found that the hard-core boson sys- 
tem with competing interactions on anisotropic triangu- 
lar lattices is locked to a series of commensurate super- 
solid phases for finite size systems, separated by series 
of jumps. This behavior is expected to occur, for exam- 
ple, in cold atom (bosons) systems on finite size optical 
lattices. In the thermodynamic limit, however, the or- 
dering wavevector Q(r/) becomes a continuous function 
of n, leading to a smooth crossover between a continuous 
set of novel incommensurate supersolids phases. 
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